A Krylov subspace algorithm for multiquadric interpolation in many dimensions

We consider the version of multiquadric interpolation where the interpolation conditions are the equations s(xi ) = fi , i = 1, 2, . . . , n, and where the interpolant has the form s(x) = ∑n j=1 λ j (‖x − x j‖ + c2)1/2 + α, x ∈ Rd , subject to the constraint ∑nj=1 λ j = 0. The points xi ∈ Rd , the right-hand sides fi , i = 1, 2, . . . , n, and the constant c are data. The equations and the constraint define the parameters λ j , j = 1, 2, . . . , n, and α. The resultant approximation s ≈ f is useful in many applications, but the calculation of the parameters by direct methods requires O(n3) operations, and n may be large. Therefore iterative procedures for this calculation have been studied at Cambridge since 1993, the main task of each iteration being the computation of s(xi ), i = 1, 2, . . . , n, for trial values of the required parameters. These procedures are based on approximations to Lagrange functions, and often they perform very well. For example, ten iterations usually provide enough accuracy in the case d = 2 and c = 0, for general positions of the data points, but the efficiency deteriorates if d and c are increased. Convergence can be guaranteed by the inclusion of a Krylov subspace technique that employs the native semi-norm of multiquadric functions. An algorithm of this kind is specified, its convergence is proved, and careful attention is given to the choice of the operator that defines the Krylov subspace, which is analogous to pre-conditioning in the conjugate gradient method. Finally, some numerical results are presented and discussed, for values of d and n from the intervals [2, 40] and [200, 10 000], respectively.